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Abstract 

We consider the product of n complex non-Hermitian, independent random matrices, each of 
size N x N with independent identically distributed Gaussian entries (Ginibre matrices). The joint 
probability distribution of the complex eigenvalues of the product matrix is found to be given by a 
determinantal point process as in the case of a single Ginibre matrix, but with a more complicated 
weight given by a Meijer G-function depending on n. Using the method of orthogonal polynomials 
we compute all eigenvalue density correlation functions exactly for finite N and fixed n. They are 
given by the determinant of the corresponding kernel which we construct explicitly. In the large- N 
limit at fixed n we first determine the microscopic correlation functions in the bulk and at the edge 
of the spectrum. After unfolding they are identical to that of the Ginibre ensemble with n = 1 and 
thus universal. In contrast the microscopic correlations we find at the origin differ for each n > 1 
and generalise the known Bessel-law in the complex plane for n = 2 to a new hypergeometric kernel 
oF n -i. 



1 Introduction 



The topic of random matrix theory (RMT) enjoys an increasing number of applications in physics, 
mathematics and other sciences, and we refer to pQ for a recent compilation. This statement holds 
both for problems with real eigenvalues as well as with complex eigenvalues. In the latter case the 
classical ensemble was introduced by Ginibre [2] who considered complex non-Hermitian matrices X 
of size N x N, with all matrix elements having independent normal distributions. However, for certain 
applications it is not sufficient to introduce a single random matrix, e.g. when considering transfer 
matrices. The problem of studying such products of random matrices is as old as RMT itself and was 
introduced by Furstenberg and Kesten [3]. Applications for products of random matrices range from 
combinatorics [4] over Lyapunov exponents in diffusion problems [5] to Quantum Chromodynamics 
[6j to name a few recent works. 

One problem one has to face is that the product matrix often has less symmetry than the individual 
matrices. For example the product P of two Hermitian matrices is in general no longer Hermitian, 
and thus acquires a complex spectrum. Therefore in the literature two types of products have been 
considered, those which are Hermitised by considering P'P with real spectra, and those which are 
studying the complex eigenvalues of P itself. In this paper we will consider the latter problem in the 
simplest setting, by multiplying a fixed number n of independent Ginibre matrices: P n = X1X2 ■ ■ ■ X n . 
This problem has found a renewed interest in recent years in the mathematics and physics literature. 
In particular the mean or macroscopic density of complex eigenvalues given by the re-th power of 
the circular law was derived independently using diagrammatic methods [7J, singular values [8] or 
empirical spectral distributions Moreover it was shown that the product of a wide class of 
non-Hermitian Wigner matrices with independent identically distributed (iid) non-Gaussian entries 
has the same spectral distribution in the large- N limit as the product of Gaussian matrices [10j . The 
same macroscopic universality was observed numerically in [9] multiplying matrices from different 
symmetry classes. 

The aim of the present paper is to provide more detailed information for the product matrix, by 
exactly solving the problem at finite matrix size N first, using orthogonal polynomials in the complex 
plane. This goal was previously achieved for n = 2 complex non-Hermitian [6], quaternion real pj-] and 
real asymmetric matrices [12] . These can be considered as non-Hermitian generalisations of Wishart- 
Laguerre ensembles, see [13] for a review. From [6] we will borrow the idea of parameterising the n 
matrices Xj in order to compute the joint probability distribution function (jpdf) of P n . Our approach 
enables us to identify various large- A limits and to resolve the fluctuations in 3 different regions: the 
origin, bulk and edge. A new origin limit was previously found in [6] for n = 2, see also [H] for gap 
probabilities in that ensemble. The limiting microscopic density at the edge for n = 2 was shown in 
[15J to coincide with that of the Ginibre ensemble n = 1 [16\ I17j. and was conjectured to universally 
hold for arbitrary n, based on numerics and including rectangular matrices with various symmetries 
as well [9]. The same universal edge scaling limit can also be found for truncations of unitary matrices 
[18] . sub-unitary matrices [19], and for the induced Ginibre ensemble [20] , 

We find that at the origin each n > 1 leads to a new class of hypergeometric kernels, whereas 
in the bulk and at the edge the correlations agree for all n with the Ginibre ensemble and are thus 
universal. Our findings add another facet to the property of universality of correlation functions in the 
complex plane. This question was answered for a single matrix, either with non-Gaussian invariant 
distributions in the bulk [21J, or with non-invariant iid elements [22j in the bulk and at the edge. 

The paper is organised as follows. In the next section [2] we briefly list the main results obtained in 
this paper. A detailed derivation of the jpdf is given in section[3l including our matrix parametrisation 
and the computation of the weight function in the complex plane. Section 0] provides the complete 
solution for finite- A at fixed n for all eigenvalue density correlation functions in terms of orthogonal 
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polynomials and their kernel. The following sections are dealing with the large-iV limit, first for 
the density at large matrix size and large argument in section [5j Section [6] is devoted to the three 
microscopic limits mentioned earlier, before concluding in section Several technical details are 
reported in the appendices IA1 - |Pl 



2 Main results 

Given the product P n of n independent matrices Xj, j = 1, . . . , n, each of size N x N drawn from the 
Ginibre ensemble with Gaussian distribution proportional to exp[— TrX^Xj], 

P n = X x X 2 ...X n . (2.1) 

The partition function Z n can be expressed as an integral of the joint probability distribution function 
Vjpdf of the complex eigenvalues Zj, j = 1, . . . ,N of P n , given by 

. N N . N 

Z n = C n U (d 2 Z aWn (z a )) \z b - Z a f = / J] d 2 Z a V Jpd f({z}) , (2.2) 
^ a=l b>a a=l 

where C n is some known constant. The weight function w n (z) that depends only on the modulus is 
given by the so-called Meijer G- function 



w n {z) = tt"- 1 ^ 



*r • (2-3) 



The corresponding kernel of polynomials orthonormal with respect to that weight reads 

_N-1 

(nk\ 



N-l 

The Appoint density correlation functions then easily follow to be the determinant of that kernel, 

R^\ Zl ,...,z k ) = /jv^'mi Y~ J d 2z k+x ■ ■ ■ d 2 z N V jp df({z}) (2.5) 

K^\zi, zj)] . (2.6) 



det 

l<i,j<k 



After briefly summarising our results for both finite n and iV which are new for n > 2 we turn to the 
large- N limits. 

For large N and large arguments \z\ 3> 1 the eigenvalue density behaves as 
From this the mean or macroscopic large-iV density can be read off: 



2.-2 

lim ±-N n R^ } (z = N n / 2 w) = 9(1 - \w\) . (2. 



ptlroH = lim ^N n R[ n \z = N n / 2 w) = ^ 



This result was previously derived in [7J [HJ [10] and we have rescaled variables such that the support 
is now the unit disc. Our new results are for the microscopic limit at the origin, in the bulk and at 
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the edge of the spectrum. We first quote the latter, by zooming into the region around the edge of 
the support, which is z ^ in eq. (|2 . T() : 

P^dgeiO = ^™ 0O " iVn " li2 i n) (z = N n /\1 + £^fi/N)e i *) = -Urfc(v/2£) • (2.9) 

This result depending only on the radial distance from the edge is universal in the sense that it agrees 
for all n with that of the Ginibre ensemble [16} I17j with n = 1 (see also |23j). and the non-Hermitian 
Wishart ensemble [15] with n = 2 derived previously. 

For the bulk scaling limit we first have to unfold, in order to have a flat density locally. In our 
case this can be obtained by a simple rescaling £ = ^/riz 1 ^ n , 

flMfc) = \f-lR^ (z = , (2.10) 



and correspondingly for higher correlation functions. Taking the scaling limit 1 <C £ < to be in 
the bulk we obtain the following answer for the asymptotic kernel 

2V-»oo;|£ij|»l 



i / ae* 

1 — ' exp 



(l-n)/2 

(2.H) 



5(l6l 2 + fe| 2 + 6^) 



This is equal to the rescaled universal kernel in the bulk of the Ginibre ensemble n = 1, in which the 
bulk and origin limits coincide. The phase factor in front of the exponential is irrelevant as it drops 
out after taking the determinant in eq. (I2.6p . 

Finally let us emphasise that there exists a third microscopic large- iV limit at the origin which 
differs for every n > 1. Its kernel is simply given by extending the sum in eq. (|2.4f) to infinity, leading 
to a hypergeometric function: 

K %rigin( z i> z j) = J im K n\ z ^ z 3) = \J w n (zi)w n (zj) Q F n -i(-; 1; ziz*) ir~ n ,n> 1 . (2.12) 

Here we have \zi\ = 0(1). For the Ginibre ensemble n = 1 the limiting kernel is given by that of 
the bulk limit above, whereas for n = 2 one has a iT-Bessel function for the weight times an 7-Bessel 
function from the sum in the kernel [6j, which was also investigated in [13] including its Fredholm 
determinant. 



3 Derivation of the joint eigenvalue distribution 
3.1 Outline of the method 

We are interested in statistical properties of eigenvalues of the product 

Pn = X1X2 ■ ■ ■ X n (3.1) 

of n independent Ginibre matrices of dimensions N x N. The partition function of independent 

matrices factorises into a product of independent partition functions for each matrix. We however 
write it in a compact way 

Z n = f[ I \DXj\ e~ TrX l Xj = f \DX\ e^^i^i . (3.2) 
j=i J J 
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because we are interested in the eigenvalues of P n . In this notation DX is the Euclidean volume form 
(external product of all independent one- forms) and \DX\ is the corresponding unoriented volume 
element: 



n n N 

\DX\ = \[\DX j \ = J{ [J ((dReX i ) a6 (dImX i ) a6 /2 7 r) 

j=l j=la,b=l 
n N 

=n n . 

j=l a , 6=1 



1,2, 



(3.3) 



, n and 



The star denotes the complex conjugate. The index j runs over the set of matrices j 
a and b over rows and columns of the matrix a = 1, . . . , N and b = 1, . . . , N. 

Our first goal is to rewrite the integrand of f)3.2[) in new variables that can be divided into two 
classes: the first one comprising variables that are directly related to the eigenvalues of P n and the 
second one comprising variables independent of the eigenvalues of P n . Having done that we can try to 
integrate out the variables of the second class to obtain an integrand depending only on eigenvalues 
of P n . We obtain in this way the joint probability distribution function Vjpdf f° r eigenvalues of P n . 
We basically follow the Dyson procedure but introduce some pertinent modifications that enable us 
to handle the case of eigenvalues of the product matrix P n in eq. (|3.ip . 

It is convenient to consider a cyclic block matrix B, following |7]: 



B 



( 


*1 





. 











x 2 


. 














A 3 ■ 














. 


• • x n - 


V x n 








. 






(3.4) 



/ 

Each block Xj in B is of dimensions N x N, so the whole matrix has dimensions nN x nN . One can 
easily see that the re-th power B n has identical eigenvalues as P n , but each being n times degenerate. 
Indeed, B n is a block diagonal matrix. The first diagonal block is X\X 2 ■ ■ ■ X n , the second one is 
X2X3 . . .X n Xi, and all others are products of consecutive cyclic permutations of X\, X2, ■ ■ ., X n . 
Each diagonal block of B n has the same eigenvalues as P n , so therefore each eigenvalue of P n appears 
in B n in n copies. 



3.2 Change of variables 

In this section we change variables, as outlined above, and find the corresponding Jacobian. Without 
loss of generality we illustrate our method for n = 3. We use the following parametrisation of the 
block matrix: 

/ X l \ 
B = 0X 2 \ 
\ x 3 / 

/ Ut \ / Ai + Ti 
= U 2 

V U 3 ) \ A 3 + T 3 

This is a generalisation of Osborn's idea who considered the chiral block matrix for n = 2 [6j . It is also 
known as QZ or generalised Schur decomposition, see e.g. |24j section 7.7. The diagonal block unitary 



\ fu^ \ 

A 2 + T 2 U 2 l 

/ V o^- 1 / 
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matrix U on the left hand side has three diagonal unitary blocks Uj of dimensions N x N restricted 
to the right coset U(N)/U(1) N . Each block Aj is a N x N diagonal matrix Aj = diag(xji, . . . ,Xjff)- 
Each block Tj is an N x N strictly upper triangular matrix. One should note that the eigenvalues of 
B 3 (and thus also of P 3 ) are equal to x\ a x 2a x 3a for a = 1, . . . , N . 

In the original representation, the block matrix B is parametrised by the elements of the matrices 
Xj. Each matrix Xj has N 2 elements, and each element is a complex variable (two real degrees of 
freedom). So altogether there are 6iV 2 real degrees of freedom. This number matches the number 
of degrees of freedom of the parametrisation on the right-hand-side: each Uj has N(N — 1) real 
parameters, each Aj has 2N and each Tj has N(N — 1). It makes 2N 2 times 3, that is 6N 2 real 
degrees of freedom as in the original parametrisation. In Appendix A we calculate the Jacobian for 
the change of variables (|3.5p and integrate out the Tj's and Uj's. The corresponding integrals factorise 
and the constants can be read off from the known Ginibre ensemble. We are left with an integral 

f " N N N 2 1 

Z3 = Cs / II II ^ X i a II \ X ^> X ^Zb - X la X 2 aX 3a \ 2 e"^ 1 £°=! lX:ial , C 3 = (^Q^.-n N)3 ' 
-> j=\a=l b>a ^ ^ ' 

N 1 (3 ' 6) 

that depends on eigenvalues of the three matrices in the product. Here G[N + 1] = Yl k=0 k\ is 
the Barnes G-function. As discussed the z a = xi a X2 a X3 a , a = 1,...,N correspond to the complex 
eigenvalues of the matrix P n . So we can now rewrite the last formula by integrating out the Xj a 's, and 
by expressing the result in terms of the z a 's. This gives the standard form of the partition function 

r N N 
Z 3 = C 3 / J] {d 2 z a w 3 (z a )) H \z b - z a \ 2 , (3.7) 

J 0=1 b>a 

however with a non-standard weight function: 

W3 (z) = [ d 2 Xl d 2 x 2 d 2 x 3 6&(z - x lX2 x 3 ) e-d^l'+l^^+l^sl 3 ) . ( 3 . 8 ) 



3.3 The weight function 

It is easy to make the generalisation of the previous two subsections to the product (13. If) of any number 
of matrices n. Apart from the constants the only difference is that the weight function w 3 (z) has to 
be replaced by w n (z): 



N N 



n-C n U (d 2 z a w n (z a )) J| \z b - z a \ , C n - cjyi qi n , i] n N)n 
o=l b<a \ • L J > 



(Ci) n , (3.9) 



that is given by: 



w n {z) 



Y[d 2 xj 5& [z-Jlxj] e-^l^l 2 . (3.10) 
i=i \ 3=1 / 

If one introduced a probabilistic normalisation to the weight function, the last equation could be 
interpreted in terms of the probability density function of a random variable obtained being a product 
of n independent Gaussian complex random variables - very close to the spirit of our original problem. 
Of course things are more complicated since in addition to the weights we have also the repulsion 
coming from the Vandermonde determinants. The second observation is that the moments of the 
weight function factorise into n independent Gaussian integrals 

J d 2 z\z\ 2k w n (z)=p_J d 2 x,\ Xj \ 2k e-^ 2 = d 2 x\x\ 2k e-^y = ^k\) n , (3.11) 
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since integrating out the delta constraint in A3. 10j) makes the a; ./-integrals mutually independent. 

One can find an integral representation of the weight function (|3,10p by integrating out one variable, 
say x n = z/{x\ . . . x n -\), and getting rid of the delta function constraint. Equivalently it can be done 
using the formula (|3.1ip for k = in polar coordinates: 



r exp|z<; 



n ■ ■ ■ r n _i ri • • • r n _i 



n > 1 , 



(3.12) 



where r,- = \xA, r = \z\. This leads to the following transformation: 



J d 2 zw n (z) 



d 2 x\ . . . d 2 x n exp 



2tt 



drr(2n) 



n-l 



drm . . . dr n _ x r n ^x 



(n ••■ r n -i) 2 



cxp 



n-l 



: y r 2 

(ri • • • r n _i) 2 ^ 3 



with 



wJz) = (2tt) 



n _ x f 00 dri... dtn-i 
Jo ri • • • r n _i 



■ exp 



(ri • • -r„_i) 2 



n-l 



n > 1 . 



(3.13) 



(3.14) 



Obviously the weight given by an (n— l)-fold integral only depends on the modulus \z\. It immediately 
follows from (|3.14j) that 



w n+ i(z) = 2ir 



dr 



w n (z/r) exp [— r 2 ] , n>0 



(3.15) 



with 1^1(2:) = e~' z which is the Ginibre weight. This recursion relation can be solved for w n +i(z) 
by using the Mellin transform in the variable R = \z\ 2 . Defining £l n (R) = w n (^/R) we calculate the 
corresponding transform: 

M n {s) = dR R s - l fl n (R) . 
Jo 

The recursion relation (13, 15ft takes the form 

n n+l (R) = 7v r^n n (R/ P )e-p , 

Jo P 



(3.16) 



and its Mellin transform factorises 



M n+l {s) = itM n {s)T{s) , 



where T(s) = f °° dpp s 1 e p . Together with the initial condition 

/>oo 

Mi(s) = dR R^expl-R] = T(s) 
Jo 

M n {s) = it n - l T n {s) . 



we thus have 
We can now write the inverse transform 



(3.17) 

(3.18) 

(3.19) 
(3.20) 



1 rc+ioo / \ 

Wn ( z ) = n n (\z\ 2 ) = n n - 1 —J^ r n (s)\z\- 2s ds = n n - 1 G% n h \z\ 2 j, (3.21 
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where the symbol Gp™(. . . \z) denotes Meijer's G-function [25], and is a string of n-zeros. Only the 
subclass Gq® of these functions, corresponding to the inverse Mellin transform of T n (s), appears in 
this context. Note that our weight functions differ from the Meijer G- functions Ggq appearing in the 
distribution of Fuss-Catalan numbers [3j. 

For n = 1 this result reproduces by construction the weight for Ginibre matrices (or by inverse 
Mellin transform 



w\(z) = exp [- 



"01 I 



(3.22) 



For n = 2 the integral eq. (|3.14p can be performed analytically, leading to a modified Bessel-function 

EI: 



w 2 (z) = 2itK {2\z\ 



11 (jr 2 \ 0,0 



(3.23) 



Alternatively one can use the connection of Meijer G- functions to special functions in [25]. For n = 3 
we either keep a single integral over the previous .fT-Bessel function, or we write the result explicitly 
in terms of the Meijer G-function: 



ws(kl) 



An 2 



dr 



K (2\z\/r)exp[ 



71 ^03 \ 0,0,0 



(3.24) 



For an alternative representation in terms of hypergeometric functions and series of polygamma special 
functions we refer to [26j section V case f). 



4 The orthogonal polynomial approach 

Since the partition function Z n (|3.9ft expressed in terms of the eigenvalues of the matrix P n has the 
standard form of the product of the weights times the absolute value square of the Vandermonde 
determinant we can use the standard technique of orthogonal polynomials in the complex plane to 
determine the fc-point correlation functions [23]. For recent reviews on non-Hermitian random matrix 
theory and its corresponding polynomials we refer to |18} 113]. We refer the interested reader to these 
references for details. 



4.1 Orthonormal polynomials and kernel 

From the fact that our weight w n {z) is angle- independent it immediately follows that the corresponding 
orthogonal polynomials are monic Pk(z) = z k . Indeed, the integral 



pill f'OG 

d 2 z w n (z)z k z* 1 = / ttye^*-') / dr rw(r) = 5 kl h^ ] (4.1) 
Jo Jo 

gives zero for k ^ I. Using eq. ()3.11j) we immediately find that the squared norms are 

= (7Tfc!)« , (4.2) 
so that the corresponding kernel of orthonormal polynomials reads 

N-l 

K%\z hZj ) = JwfrMzj) £ ^(^f • (4.3) 

k=o {7rK ' ) 

In view of the results later in section [6751 we call this a truncated hypergeometric kernel, just as in the 
Ginibre case n = 1 this is a truncated exponential. 
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Following [23] this kernel determines the fc-point correlation functions: 



R { k\zi,...,z k ) = det 

i<j j,<fc 



(4.4) 



In particular, the one-point correlation function or spectral density is R\ 

point correlation function is R^\z,u) = K^\z, z)K^\u,u) - h 
going to discuss these in the large- N limit in the next section. 



(*), 



K 



K^ 1 ' (z,u)K^\u, z) etc. We are now 



(n). 

AT 1 



the two- 



5 Eigenvalue density at large- N 

In this section we study the eigenvalue density for large- A (and large argument). The purpose of that 
is two-fold, as we will not only find the mean or macroscopic density after a suitable rescaling. The 
expression for large but finite- A that we will derive below enables us to identify the different regions 
in the complex plane where we can take different microscopic limits that describe the fluctuations in 
that region. 

We begin with the definition of the spectral density or one-point correlation function in terms of 
the kernel: 

N-l | |2fc 

R^\z) = K%\z,z) = w n (z) J2fl^ = M*)T n {\z\\N). (5.1) 

k=o [ K -> 

We note that according to its definition eq. (|2.5p it is normalised to the number of eigenvalues, 
f d 2 zR^\z) = N, rather than to unity. 

Now we are going to determine the behaviour of R^ (z) for large N and for large \z\. The 
asymptotic large- behaviour of the weight w n (z) can be taken from [27J. It is rederived here for 
completeness in Appendix |B] (see equation (IB. 5ft ). The large- N behaviour of the truncated sum 
T n (x,N) for x = \z\ 2 of order N n was already studied in [28], and we give some details of the 
derivation in Appendix O (see equation (1C.6P ). Putting together (|B.5P and (1C.6P we obtain the 
following behaviour for large- N and for \z\ < N n ^ 2 : 



R ( i\z) - ^rfc [ v '^, M I . (5.2) 



M"" 2 1 r f Vn{\z\ 2/n -N) 



rnr 



2 I y/2\z\V n 



It is convenient to recast this result into a rescaled density with compact support that is normalised 
to unity. Using the rescaled variable w = zN~ n / 2 we define the following density, for which the radius 
of the eigenvalue support approaches one for TV" — > oo: 



——2 

lim N n ~ l Rf\N n / 2 w) = -erfc 



(w = 




2 I V2\w\y n 



(5.3) 



The factor A^" in the first step comes from the rescaling of the volume element d 2 z = N n d 2 w, and the 
additional factor A -1 from the normalisation to unity. In the second step we have simplified p n (w) 
by Taylor expanding the argument around unity, \w\ = (2/n)(\w\ — 1) + 0((|w| — l) 2 ). This is because 
the complementary error function changes only in a narrow strip around the unit circle \w\ = 1, of a 
width proportional to 1/ y/N. So far we have derived a finite size form of the eigenvalue distribution 
of the product of n independent Ginibre matrices. This form was known for n = 1 [16} [T7] and n = 2 
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[15] . It was conjectured to hold for any n and tested numerically in [9] except that the dependence of 
the coefficient inside the erfc-function on n was unknown. We have filled this gap here by deriving eq. 
(I5.3j) for any n. We see that the width of the crossover region around the edge |w| = 1 is proportional 
to the square root of the number of multiplied matrices \fn. 

It is instructive to compare this result with the limiting density for large-A of the n-th power 
of a single Ginibre matrix. This density is derived in Appendix [D] and is given by exactly the same 
distribution, however with a different dependence on n: 



Pn(w) = J-! -erfc (\w\ - 1) . (5.4) 




Here, the width of the crossover region is proportional to n and not to y/n. In a sense, the finite size 
corrections are stronger for the n-th power than for the product of n independent Ginibre matrices. 

Actually the same seems to be true for a broader class of matrices called isotropic [29]. Isotropic 
matrices are known to have the property that the product of iid matrices has the same limiting 
distribution as the corresponding power of a single one, as was discussed in [29]. For this class the 
finite size formula is not known but numerically one also observes that the product of independent 
matrices converges faster to the limiting distribution than does the corresponding power [29J. 

We close this section by rederiving the mean or macroscopic density in the limit N — > oo. It is 
obtained by taking the limit N — > oo of (|5.3|) : 



I I--2 

PmLro( W ) = J im Pn(w) = 6(1 - \w\) , (5.5) 



where 0(x) is the Heaviside step function. This result was derived in [7] using planar diagrams and 
rigorously in [8] [10] . We also refer to [8] for an extensive list of references regarding the mathematical 
literature on products of random matrices. 



6 Microscopic large- N limits and universality 
6.1 Universal edge limit 

In this subsection we analyse the fluctuations of the eigenvalues around the edge of the spectrum. We 
will only consider the density itself, higher order correlation functions can be dealt with analogously. 
Our starting point is the density at large-A and large-|z|, eq. (|5.2p . As was already mentioned in the 
previous section the complementary error- function rapidly vanishes for \z\ > N n ^ 2 , on a strip of width 
1/y/N. We therefore introduce a new scaling variable £: 

z = N n ' 2 (l + ? v 7 ^) e # • (6.1) 

The additional rescaling with \Jn is made in order to normalise properly. We can therefore define the 
microscopic density at the edge by 

lim P%L(S) = lim nN n - l R^ ] (z = A"/ 2 (l + iy/^/N)^) = -Urfc(V2~0 , (6.2) 

where for each real dimension the density R± has to be rescaled by y/nN n l 2 . The resulting expression 
agrees with the n = 1 result for the Ginibre ensemble derived in |16} [T7] (see also [23J ) . For n = 2 
the same result was previously derived in [15]. We have shown here that the result holds true for any 
fixed n and is thus universal. This behaviour was already observed numerically in [9]. 
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6.2 Universal bulk limit 

In this subsection we discuss the behaviour of the /c-point correlation functions J?£ (21, . . . , z k ) defined 
in (|4.4|) . in the limit N — > oo for the following scale separation 1 -C \z k \ < N n / 2 . By comparing to eq. 
(|5.2p this separation will guarantee that we are neither close to the origin to be discussed later, nor 
to the edge limit discussed above. 

In principle the microscopic bulk limit is obtained in three steps. We would first have to rescale 
z = N n l 2 w in order to get a compact support. In the second step we would zoom into the fluctuations 
which are in this case of the order l/iV n / 2 , w = v/N n / 2 . In the last step we would unfold the variables 
in order to get a locally flat density, with a mean level spacing of unity. Because of the macroscopic 
density being given by 

2 -G(l - M) , (6.3) 



2 



An) / x 



nir 

the unfolding would read £ = yjnv l l n . Because the rescalings from the first two steps compensate 
each other here, we can do all steps in one. Our microscopic scaling limit is thus given by the change 
of variables 

K |£ = V^ 1/n | < VN , (6.4) 

which will become ^-independent in the large- N limit. Note that in the variable z 1 ^ the fluctuations 
are again of order 1/y/N as in the Ginibre ensemble. The change of variables for the /c-point densities 
can be best seen by considering the two-dimensional differential 

d 2 i = n- 1 \z\l- 2 d 2 z . (6.5) 

We thus obtain for the microscopic density 

R^(0 = n\z\ 2 -T,R^ ( z = {aV^T) = n\z\ 2 -l±-\z\l- 2 = - , (6.6) 



which is constant on a disk of radius yN and zero otherwise. So indeed the eigenvalues are uniformly 
distributed in the £-plane for large |£| inside the support. More generally, the microscopic /c-point 
correlation functions are defined as 

= ■ ■ ■ z k \ 2 - 2/n R^\z 1 , ...,z k ), (6.7) 

where Zj = (£,j/\/n) n on the right hand side, from eq. (|6.4|) . All /c-point correlation functions 
R k n \^i, . . . , can be expressed in terms of a redefined kernel 



4 n) (&>--->&)= det Kfi'fatj) , (6.8) 
that is given by 

=n\z i z j \ 1 - 1 ' n K$\z i ,z j ) = n 2 - n \^r l K^ {fa/y/n) n , ) . (6.9) 
In the limit N — > oo and 1 <C < the limiting kernel in the bulk is given by 

^^^^^-^^^^/^^^(^/^^^((^/n)") , (6.10) 
where we have for asymptotically large argument (|B.5|) 

lim ((fA/n)") = 4= (2vr 3 n) ( ^ 1)/2 ^""eH*!' , (6.11) 



-(")/ 
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and for asymptotically large argument and large- N (jC.6h 



lim T n ((£i£/n) n ) 



7T\/ TL 



(2vr 3 /n) 



(l-n)/2 



(6.12) 



Putting all together we obtain 



i / 



* \ (l-n)/2 



exp 



(6.13) 



We see that the kernel is Hermitian if £?)* = K^^j,^). It is equal to the kernel of the Ginibre 
ensemble multiplied by powers of the phase factors e 1 ^ = The kernel (|6,13p is unitarily 

equivalent to the kernel of the Ginibre ensemble. When one calculates the correlation functions 
(&!••■ >£fc) all phase factors cancel and one obtains exactly the same /c-point correlation 

functions as for the Ginibre ensemble with n = 1. Therefore all correlation functions of the product 
matrix P n in the bulk limit are universal. In particular the two-point correlation function is given by 



R 



l 



71" 



1 _ e -l£i-6l : 



3.14) 



which describes correlations of eigenvalues at distances |£i — £2! of order unity. 



6.3 Microscopic origin limit: n classes 

In this subsection we investigate the remaining region where a non-trivial microscopic limit can be 
obtained, that is the vicinity of the origin. While the microscopic bulk and edge limits yield universal 
results for arbitrary fixed n, at the origin we obtain a different kernel for each n. For n = 1 it coincides 
with the bulk limit whereas for n = 2 we obtain a kernel containing Bessel-iT and Bessel-J functions 
which was known based on the work [6], see also |14| for more details. 

If we look back at the previous two limits the edge correlations were found for \z\ ~ N n l 2 or \w\ ~ 1 
in rescaled variables z = N n / 2 w, with fluctuations of order l/y/N. The bulk limit was obtained by 
considering 1 <C \z\ < N n l 2 or < \uu\ < 1, with fluctuations again of order 1/y/N in terms of the 
rescaled variable £ = y/nz l l n . Here we will take the limit N — > 00 keeping \z\ = 0(1), or \w\ ~ N~ n / 2 . 
Because of that the scaling is very simple in z: the weight w n (z) which is iV-independent remains 
unchanged, and the sum inside the kernel eq. (15. ip simply has to be extended to infinity: 

K wiain(. Z i> Z j) = lim K N\ Z hZj) = \/w n (zi)w n (zj) F n -i (-; 1 ; ZiZ* ) 7r~ n , n > 1 . (6.15) 
a N—too v 

Here we have used that the infinite sum eq. (|C.7P is related to a generalised hypergeometric function, 

T ^ x ) = E T^yT = ^ Q^-iH ?! x) , (6.16) 
for n > 1. The vector 1 contains n — 1 elements of unity. For n = 1 we simply have 

r i( x ) = E^ = - ex PW • (6-17) 

fc=0 

It cancels the exponent exp[— \z\ 2 ] coming from the weight and thus leads to a constant density. 
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For n = 2 the result can also be expressed in terms of elementary functions, due to 



oo 



k 




E 



x' 



(6.18) 



fe=0 



(ttA:!) 2 



7T 



This leads to the known kernel (see e.g 



[14J) 




(6.19) 



in terms of Bessel functions, with the density P or igi n (z) — 

--Ko(2|2:|)-fo(2|,z|). For n > 3 our ker- 
nel eq. (|6.15p seems to be new. Although other random matrix ensembles with kernels containing 
hypergeometric functions exist, we are not aware of any example containing o-^n-l- 

7 Conclusions and outlook 

In this paper we have investigated the eigenvalue correlations for the product of n independent non- 
Hermitian random matrices. Explicit results were given for any finite matrix size N for all density 
correlation functions in terms of the determinant of a kernel of orthogonal polynomials in the complex 
plane. The corresponding weight functions were found to be so-called Meijer G-functions depending 
on n, whereas the joint probability distribution function remained the absolute value squared of the 
Vandermonde determinant which is standard for this symmetry class. 

So far such results were only known for up to n = 2 (in fact for rectangular matrices with an 
elliptic law). In the large- N limit we investigated three different microscopic limits where the local 
fluctuations were zoomed into: the microscopic edge, bulk and origin limit. In the former two cases 
we found a complete agreement with the Ginibre ensemble of a single matrix n = 1, after unfolding. 
The edge and bulk limit are thus universal, as it was conjectured earlier for the behaviour at the 
edge depending on the complementary error function. In contrast at the origin each new matrix in 
the product adds a new class of correlation functions, given by an n-dependent hypergeometric kernel 
together with the Meijer G- function from the weight. 

Several generalisations of the results presented here are conceivable. First of all it should be pos- 
sible to consider the product of rectangular matrices instead. We expect that the results for the edge 
and the bulk remain unchanged, whereas the origin limit will be generalised as it is already known for 
n = 2. Matrices from other symmetry classes could be considered, such as real quaternionic non-self 
dual or real asymmetric matrices. Again we expect the edge and bulk behaviour to be unchanged, with 
local microscopic modifications along the real and imaginary axis as well as at the origin. Finally also 
for each matrix more general distributions than the Gaussian Ginibre distribution could be feasible. 
In particular it would be interesting to know if the microscopic properties we computed continue to 
hold for the product of non-Hermitian Wigner matrices with iid but non-Gaussian entries. Several of 
these projects are currently under way. 



Acknowledgments: We would like to thank Eugene Strahov for useful discussions. We acknowl- 
edge partial support by the Polish Ministry of Science Grant No. N N202 229137 (2009-2012) and 
by the Grant DEC-2011/02/A/ST1/00119 of the National Centre of Science (Z.B.), as well as by the 
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A Computation of the Jacobian 



We now calculate the Jacobi matrix for the change of variables from X to (U,T,A). In order to 
determine the relation between the infinitesimal elements (one-forms) we differentiate both sides of 
(1331). We obtain 



dXi 
<1X (0 dX 2 






U2 1 
[/3- 1 



(A.l) 



where 
and 



dYj = dAj + dTj + dMj , (A.2) 



(dMi) ab = {dA x ) ab zit - zia (dA 2 ) ab + £ {dA x ) ac {T x ) cb - J2 {Tx) ad {dA 2 

c<b d>a 



(dM 2 ) ab = (dA 2 ) ab x 2b - x 2a {dA 3 ) ab + £ (dA 2 ) ac (T 2 ) cb - £ {T 2 ) ad {dA 3 ) db , {A 3) 

c<b d>a 

(dM 3 ) ab = (dA 3 ) ab x 3b - x 3a (dA^ + J2 ( dA z)ac (Ts) cb - Yl ( T ^ad (dAJa . 

c<b d>a 

The dAj = UJ x dUj are anti-Hermitian matrices with zeros on the diagonal. One can choose {dAj) b 
and (dAj)* ab in the upper triangle (a > b) to be independent infinitesimal elements of dAj (or more 
precisely - independent one-forms). The elements (dAj) ab in the lower triangle (a < b), can be 
expressed by those in the upper one: (dAj) ab = — (dAj) ba . The diagonal elements of dAj are zero 
since Uj are restricted to U(N)/U(1) N . 

Now we are ready to calculate the Jacobian. Let us do it gradually. First observe that the Jacobian 
for the change of variables dX = \JdYU~ 1 eq. (jA.ip equals one since the matrix U is unitary. So we 
have \DX\ = \DY\. Now we can calculate the Jacobian for the change from dY to (dA,dT,dM) in 
eq. (|A.2p . To this end let us write down (|A.2|) in an explicit index notation 

{dYj) gjb = d{Aj)ab + {dMj) ab for a = b , 

{dYj) ab = d(Tj) ab + {dMj) ab fora<6, (A.4) 

(dYj) ab = {dMj) ab for a > b . 

There are also corresponding equations for the complex conjugates. The first line, for a = b, could be 
alternatively written as {dYj) aa = dxj a + {dMj) aa since d{Aj) ab = 5 a bdxj a , or in shorthand notation 
dA = dx. As one can see from (|A.4p the Jacobi matrix for the linear transformation from the (iY-basis 
to the {dx, dT, (iM)-basis is upper triangular and has a Jacobian equal one, so we have 

\DX\ = \DY\ = \Dx\\DT\\DM\ . (A.5) 

The measure DM in the last equation is a product of one-forms {dMj) ab for a > b. They can be 
expressed as a linear combination of the one-forms dA's in eq. (|A.3P : 



\DX\ = \Dx\\DT\\DA\ 



dM 
~dA 



(A.6) 
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Let us calculate the Jacobian \J\ = \dM/dA\ for this transformation (|A.3j) . It has a specific algebraic 
structure that imposes a certain ordering of the indices that appear in (|A.3p . In fact this ordering 
implies that the Jacobi matrix can be written as an upper triangular block matrix. To see this let 
us first recall that we are considering only the sector a > b corresponding to the lowest line of (|A.4p . 
Second of all, let us note that as a result of the Tj's being strictly upper triangular matrices we have 
further inequalities for the indices in ()A.3p : 

a> b> c and d> a> b . (A. 7) 

An inspection of the indices of the one-forms appearing in (|A.3[) . (dMj) ab (a > b), (dAj) ac (a > c), 
(dAj)db (d > b), shows that they are all indexed by ordered pairs of indices whose first index is larger 
than the second one. It is convenient to introduce an increasing ordering in the set of such ordered 
pairs. A pair ab is said to be less than a pair cd (and denoted by ab < cd) if a > c or if a = c and b < d. 
This choice of ordering does not look very intuitive but it is convenient for our purposes. For example 
for N = 4 there are six ordered pairs and they are ordered as follows 41 < 42 < 43 < 31 < 32 < 21. 
We can now introduce a single index a = 1, ... ,6 that preserves the increasing ordering of the pairs: 
41 -> 1, 42 -> 2, 42 3, 31 -> 4, 31 -)• 5, 21 -> 6. Of course, we can do this for any N. Using the 
index a we can concisely write equations (IA.3j) as 

(dMj) a = J2(J jj/ ) aal (dA k ) al . (A.8) 
j'a' 

Now we want to argue that the Jacobi matrix is upper block triangular, that is (Jjj<) aa ' = for all 
a > a' . Clearly the matrix has diagonal blocks a = a' (IA.3j) . It also has off-diagonal blocks coming 
from the sums on the right-hand side of (|A.3|) . We now show that the terms in the sums contribute 
to the upper triangle a < a'. Indeed, as follows from the inequalities (1A.7I) the pairs of indices ab of 
{dMj)ahS on the left hand side of (|A.3|) are smaller (in the sense defined above: ab < ac and ab < db) 
than the corresponding pairs ac and db of (dAj) ac and {dAj)db i n the sums on the right hand side 
of (|A.3|) . This is equivalent to saying that the sums run over a"s such that a < a'. There are no 
terms for a > a' on the right hand side of (1A.3|) and thus the Jacobi matrix (|A.8|) (Jjj')aa' = for 
a > a'. Since the Jacobi matrix (Jjj>) aa i is upper block triangular, its Jacobi determinant |J| is 
identical to the determinant of the corresponding block diagonal matrix (Jjj')aa' that has the same 
diagonal blocks: \J\ = \J\. In this way we have reduced the calculation of the determinant \8M/dA\ 
to the calculation of the determinant \8M/dA\ of the block diagonal matrix J corresponding to the 
transformation obtained from (IA.3j) by skipping upper-triangular terms, that is the sums: 

dMi) = (dAx) ah x lb - xia (dA 2 ) ab , 

/ ab 



dM 2 \ = {dA 2 ) ab x 2b - x 2a (dA 3 ) ab , (A.9) 

/ ab 

( d ^ 3 )ab = ^ ~ ' 

The Jacobian can be calculated by the use of the external (wedge) product of one-forms (dMj) ab . Let 
us first calculate it for a given pair ab. Using (1A.3|) we have: 

(dM 1 ) ab A(dM 2 ) ab A{dM 3 ) ab = 

{dAi) ab A (dA 2 ) ab A (dA 3 ) ab (xi b x 2b x 3b - xi a x 2a x 3a ) . 
Taking also into account the corresponding independent equations for the complex conjugates, we 
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obtain a transformation of the 6-form in the whole ab sector : 

(dmU A {dM 2 ) ab A (dM 3 ) ab A (dM 3 f ab A (dM 2 )* ab A {dM{)* ab = 

{dA x ) ab A {dA 2 )ab A (dA 3 ) ab A (dA 3 )* ab A (cL4 2 )* 6 A (cL4i)* 6 (A.ll) 

I |2 

We repeat this calculation independently for all sectors a > b and eventually obtain the volume form: 

N 

\DM\ = \J(x)\\DA\ = Yl \x lb x 2b x 3b - x la x 2a x 3a \ 2 \DA\ , (A.12) 

a>b 

and 

\DX\ = \J(x)\\Dx\\DT\\DA\ . (A.13) 

We are ready to write down the partition function (|3.2|) in new variables, for general n. The 
Gaussian weight function (|3,2p assumes a form: 

e - Tv T,"=l X j X i — e -Y^,"=lJ2a = l\ X ja\ 2 e -J2"=lJ2a>b\( T j)ab\ 2 ^ (A. 14) 

that is independent of dA^s and thus of dU^s: dAj = U~ l dUj. So we can integrate out the dU- 
variables. This integration gives a constant equal to the volume of the coset U(N)/U(1) N . Also the 
integration over the (Tj) ab s can be done since it is an independent Gaussian integral. By comparing 
to Ginibre we obtain for each j 

J \DUj\ J \ DTj \ exp[-Tr T^Tj] = 1/(N\ G[N + IJtt^) = C x . (A.15) 

Denoting by C n the total constant factor coming from the integration over all dUj J s and dT^s, we 
have C n = {Ci) n . What remains is an integral over the cfe's as was claimed in section 

„ n N N 

Z n = C n / Yl II d2x i a II \ x ^X2bX 3b - x la x 2a x 3a \ 2 e" lxjal2 . (A.16) 

J j=la=l a>b 

B Asymptotic of the weight function 

In order to make our paper self contained we compute here the leading order asymptotic behaviour 
of the weights w n (z) for \z\ — > oo using the saddle point method. The result can also be found in 
Theorem 2 in reference [27] . 

We use the multidimensional representation of the weights (I3.14p : 

/n— 1 
Hidrj/rj) e~ s , n> 1 , (B.l) 
i=i 

where 

12 n_1 



The saddle point equation 



16 



has a symmetric solution n = r% = ■ ■ -r n -i = r* = |z| 1//n . The equation has no other solutions, so 
the symmetric solution is not only a local minimum of S but also the global one. At the minimum the 
function S assumes the value S* = n\z\ 2 / n . The Hessian Ha = a 92 ^ = ( i + 2S ^)\ z \ _|_ takes the 
following values at the minimum r*: H^j = 4 for i ^ j and H*ij = 8 for i = j. The determinant of 
the Hessian is det[iJ*] = n4 n_1 . So we can write an explicit formula for the leading order asymptotic 
behaviour of w n (z) for large z as 

that gives 

Wn{z) „ J_ ( 2 ^) {n - 1)/2 | z |(i-")/™ e -^l 2/n , (B.5) 



which agrees with [27] where also the error terms are computed. In particular for n = 2, 3 we have 

~ nVh-^e-W , wjW ~ ^|z|- 2 / 3 e- 3 l 2 l 2/3 , (B.6) 

v3 

where the first result n = 2 agrees with the known asymptotic of Xo(2|z|) from eq. (I3.23P [25J. In fact 
eq. (jB.5|) also holds for n = 1 where w n {z) ~ exp[— \z\ 2 ] is exact. 

The asymptotic formula (|B.5p is consistent with the recurrence relation ()3.15p in the sense that 
when one inserts the asymptotic form of w n {z) eq. (|B.5P into the integral on the right hand side of 
the recurrence relation (13. 151) one obtains the asymptotic form of w n +i(z) as given in (jB.5[) . 

C Asymptotic of the hypergeometric kernel 

In this appendix we discuss the asymptotic behaviour of the truncated sum in the kernel eq. (|4.3p . 
The asymptotic of the corresponding infinite sum was discussed in reference [28], and we again give 
some details here to be self-contained. Let us repeat the definition 

N-l k 
k=0 v ' ; 

which we would like to investigate for large N and for x of order N n , due to x = \z\ 2 < N n . We 
again apply the saddle point method. By first using Stirling's formula k\ ~ ^/2nk(k/e) k and then 
approximating the sum by an integral we obtain: 



dkk~ n / 2 e s{k) , (C.2) 



where s(k) = k \n[x] — nk \n[k] + nk. We are interested in the behaviour of T n (x, N) for large N and x. 
In this case the exponent e s ^- ) can be approximated by a Gaussian function with a maximum located 
at 

h = x l ' n , (C.3) 

being a solution of the saddle point equation 

s'{h) = ln[x] - n ln[h] = . (C.4) 
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We also have s(k*) = nk*, s"(k*) = —n/k*, so the Gaussian approximation gives: 



N 



T n (x,N)^(27T 3 )- n / 2 e nk " dkk -n/2 e -n{k-k,) /(2k«) ^ 



with k* = x l l n (|C.3|h Let us change the integration variable in the last integral to t = y / n/2k* (k — k*). 
Written in this new variable the Gaussian part assumes the form e _t , which means that values of t that 
contribute to the integral are of order one. In this narrow range of t, the factor k~ n l 2 can be treated 
as constant. Indeed, for t of order of unity and k* of order N, we have k = + t/(nk*/2) 1 ' 2 ~ k*, 
for t < (7iAs»/2)Va. The upper integration limit for t is: \fn/2k*{N — fc*). The lower one is far below 
the range of the Gaussian integrand, so it can be set to — oo. Replacing k* by x x l n (|C.3|) we eventually 
obtain 



lim T n (x,N) « — ^= (2^ 3 ) (1 n)/2 x (1 " n)/2n exp nx 1,n ^erfc I ^' 1 . (C.6) 

In the limit N — > oo the truncated sum becomes an infinite series T n (x) = ]im^^. 00 T n ^x,N). In this 
case the saddle point is always located within the range of integration, so we have 



1 / \fii{x l l n - N) 



D Eigenvalue distribution of n powers of a single Ginibre matrix 

In this appendix we discuss the eigenvalue distribution of the n-th power of a single Ginibre matrix. 
Consider a random matrix A with a spherically symmetric eigenvalue distribution p(w) = p(\w\ = r). 
Denote by p n {w) the corresponding spectral density of A n - the n-th power A, which is also angle- 
independent. If A is an eigenvalue of A, then A = A n is an eigenvalue of A n . Changing the integration 
variable in the formula below to s = r n , we have: 

f°° I' 00 27rs 2 /"~ 1 / \ f°° 
J 2wrp(r)dr = J p \^s 1/n ) ds = J 2irsp n (s)ds, (D.l) 



2—2 
rn i i / 



and thus 

Pn(r) = —p{r" n ) . (D.2) 

In particular, for the limiting macroscopic density of a single Ginibre matrix X with n = 1 we have 
Pmacro{w) = — \w\) and therefore the corresponding density of its n-power is 

i 1--2 
\w\ n 

Pmacro{w) = 6(1 - \w\) . (D.3) 

nn 

This is exactly as for the product of n- independent matrices ()5.5j) . We can also repeat calculations 
for the finite size distribution for a single Ginibre matrix in eq. (|5.3p : 

p{w) = N- l R ( l\m W ) = — erfc (V2N (\w\ - 1)) . (D.4) 

2ir V / 

Applying (|D,2p we find the corresponding distribution of eigenvalues of its n-th power: 

I ——2 

TO f rfC (^( Hl/n " 1 )) • (D ' 5) 
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It can be further simplified using a Taylor expansion at \ w\ = r = 1 to approximate the function inside 
the erfc-function: r l l n — 1 = (l/n)(r — 1) + . . . and neglecting higher order terms, as in section [5j 

\w\ n~ 2 1 / y/2N \ 

PnH = — o erfc — H-l) • (D-6) 

7rn 2 \ n / 

An important difference between this case and the power of n independent matrices is that the de- 
pendence on n in the denominator is n, and not y/n as in eq. (|5.3p . 
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